TUTORIAL 2 : NPT Lennard-Jones fluid

Authors:Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk), James Grant (r.j.grant{at}bath.ac.uk), and John Purton (john.purton{at}stfc.ac.uk)

Introduction

In the NpT ensemble number of particles, pressure and temperature (N,p,T) are kept constant, which implies that apart from the particle moves the volume is also allowed to vary. That is, an additional MC move, uniformly sampling the volume, is attempted with the acceptance probability:

\[P_{\mathrm{acc}}([\mathbf{r}_{1},V_1] \rightarrow [\mathbf{r}_2,V_2]) = \min(1, \exp \{- \beta [U(\mathbf{r}_2) - U(\mathbf{r}_1) + p_{ext}(V_{2}-V_{1}) - N \beta^{-1} \ln(V_{2} / V_{1}) ] \} )\]

where \(p_{ext}\) is the external pressure. The particle coordinates are assumed to be scaled accordingly.

NOTE: The last term under the exponent needs to be modified in the cases where, for example, linear dimensions of the (cubic) cell or \(\ln(V)\) are directly sampled instead of volume (keywords linear and log respectively).

In order to proceed to NpT simulation, we will modify a copy of CONTROL file that we used in Tutorial 1 for NVT simulation, wherein we need to introduce the directives specifying the volume move and the accumulation of additional statistics.

Navigate to the subfolder tutorial_2 (in the folder exercises). Our initial set-up contains duplicates of the input files for Tutorial 1. The CONFIG and FIELD files remain unchanged, so the initial configuration and interactions between particles will be identical to the NVT case. For your reference, the subdirectory npt-dry contains an example set of input and output files.

NPT simulation of Lennard-Jones fluid
use ortho
finish
seeds 12 34 56 78               # Seed RNG seeds explicitly to the default
nbrlist auto                    # Use a neighbour list to speed up
                                # energy calculations
maxnonbondnbrs 512              # Maximum number of neighbours in neighbour list
temperature     1.4283461511745 # Corresponds to T*=1.1876;
                                # T(in K) = T* / BOLTZMAN
pressure        0.0179123655568
steps          110000           # Number of moves to perform in simulation
equilibration   10000           # Equilibration period: statistics
                                #are gathered after this period
print           10000           # Print statistics every 'print' moves
stack           10000           # Size of blocks for block averaging to obtain statistics
sample coord    10000           # How often to print configurations to ARCHIVE.000
yamldata        1000            # collect YAML stats every 1000 move

revconformat  dlmonte           # REVCON file is in DL_MONTE CONFIG format
archiveformat dlpoly4           # ARCHIVE.000/HISTORY.000/TRAJECTORY.000 format
                                # In this case: HISTORY.000 in DLPOLY4 style
move atom 1 512                 # Move atoms with a weight of 512
LJ core
move volume cubic linear 1      # Move volume, box is cubic,
                                # linear scaling with a weight of 1
start

The directive to invoke volume moves is move volume. The parameters to specify in the present tutorial are cubic, linear and 1 for the weight of volume moves among other MC steps (see the manual for further details). Note that the weight of atom moves has been set to 512 (see move atom directive). Volume moves are more computationally intensive than single atom moves, so the rule of thumb is to attempt one volume move each time every atom has been attempted to move (so called “pass” or “sweep” through the system).

You will need to add the lines:

sample volume 1 1.0             # sample volume every V-step, with the bin size of 1 A^3

move volume cubic linear 1      # Move volume, box is cubic,
                                # linear sampling with a weight of 1

and:

pressure     0.0179123655568

Exercise

As before, in order to extract the volume sequence along the MC time from YAMLDATA.000 and plot it in gnuplot, use the commands:

[tutorial_2]$ strip_yaml.sh volume

[tutorial_2]$ gnuplot
gnuplot> plot './volume.dat' u 1:2 w l t "Volume(MC step)"

Questions to ask yourself:

• Would it help equilibration in the NpT ensemle to start the simulation after NVT equilibration?
• What would be the effect of changing pressure and/or temperature?
• How to restart the simulation for continuation?

Try to increase the pressure and run a longer simulation starting with the last saved configuration (i.e. REVCON.000).

Further exercises

Now try the extension exercises to learn more about functionality within DL_MONTE and to optimise your calculation:

TUTORIAL 2 : Exercise 1

Or move on to TUTORIAL 3 : RDF and learn how to obtain radial distribution functions (RDF).